Harmonic oscillator under Levy noise: Unexpected properties in the phase space 



Igor M. Sokolo\__ and Werner EbellingQ 
Institut fiir Physik, Humboldt-Umversitdt zu Berlin, Newtonstrasse 15, D-12489 Berlin, Germany 

Bartlomiej Dybie<H 

Marian Smoluchowski Institute of Physics, and Mark Kac Center for Complex Systems Research, 
Jagellonian University, ul. Reymonta 4, 30-059 Krakow, Poland 
(Dated: February 23, 2011) 

A harmonic oscillator under influence of the noise is a basic model of various physical phenomena. 
Under Gaussian white noise the position and velocity of the oscillator are independent random 
variables which are distributed according to the bivariate Gaussian distribution with elliptic level 
lines. The distribution of phase is homogeneous. None of these properties hold in the general Levy 
case. Thus, the level lines of the joint probability density are not elliptic. The coordinate and the 
velocity of the oscillator are strongly dependent, and this dependence is quantified by introducing 
the corresponding parameter ("width deficit"). The distribution of the phase is inhomogeneous and 
highly nontrivial. 
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I. INTRODUCTION 

Many situations in natural sciences can be successfully 
investigated adopting the stochastic level of description, 
considering the system of interest as a dynamical sys- 
tem responding to external perturbations represented by 
a noise. In the simplest situations this noise is assumed 
to be white and Gaussian. The white type of the noise is 
the consequence of the large number of independent in- 
teractions bounded in time. Its Gaussian character arises 
due to the assumption that the interactions are bounded 
in their strength. In many cases far from equilibrium the 
second assumption fails; the noise still can be considered 
as white but is now described by heavy-tailed distribu- 
tions, often of the cv-stable Levy type Such heavy- 
tailed fluctuations are abundant in turbulent fluid flows 
@, magnetizedplasmas [H, optical lattices [||, heart- 
beat dynamics [ff , neural networks J6[ , search on a fold- 
ing polymers Q , animal movement [J, climate dynamics 
Q , financial time series , and even in spreading of dis- 
eases and dispersal of banknotes ll| . These observations 
indicate the need for closer examination of properties of 
stochastic systems under such heavy-tailed noises. 

A damped harmonic oscillator under influence of exter- 
nal noise is a basic model of non-equilibrium statistical 
physics. The evolution of the state variable x is described 
by the following equation of motion 



x(t) = ~jx(t) - kx(t) + £(t) 



(1) 



(the mass of the oscillating particle is taken to be unity). 
Since the response function of the harmonic oscillator is 
known, the coordinate x(t) and the velocity v(t) = x(t) 
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are obtained by the action of a known linear operators 
on the noise term £(i), and the characteristic function 
of the joint probability distribution of (x, v) can be ob- 
tained explicitly [l2[ by generalizing the results by Doob 
[l3j. However, although the formal solution of the prob- 
lem is known for almost 30 years, it doesn't seem that 
the properties of such a basic system have attracted any 
attention of physicists or mathematicians, although they 
are astonishingly different from the ones under Gaussian 
noise. Here, we focus on the properties of the correspond- 
ing stationary distribution, which is achieved at times 
t > j/k. 

In the Gaussian case the distribution of variables (x, v) 
is a bivariate normal of the type 



p(x,v) 



1 



2tt(j x (t v 



exp 



x v 



(2) 



where o~ x and o~ v are the widths of the distributions of 
the coordinate and velocity [l_3|. Since this distribution 
factorizes into a product of x- and ^-distributions, the 
phase variables are independent. The independence of x 
and v carries over to the equipartition theorem of equi- 
librium statistical physics. Since the argument of the 
exponential is a quadratic form of the coordinate and 
velocity, the contours of equal probability density are el- 
liptic (the distribution is elliptically contoured, or simply 
elliptic). Introducing the rcscaled variables x — x/o~ x and 
v = v/a v , one defines the phase angle <fi = arctan({i/x) 
which is uniformly distributed over [— it/2, 7r/2). 

In the case of Levy noises other than the Gaussian 
one none of these properties holds. Thus, the variables 
x and v are dependent, and the dependence is stronger 
in the overdamped case (\f~h~ < 7/2) than in the under- 
damped one (\fk > 7/2). The distribution of (x, v) is not 
elliptic (see Fig. [TJ, and the phase angle shows an inho- 
mogeneous, highly nontrivial distribution (see Fig. [3}. In 
what follows we present the results of numerical simula- 
tions of the harmonic oscillator under Levy noise together 
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with analytical results corroborating these findings. We 
moreover discuss in some detail the measures which can 
be used to characterize dependence of the nonelliptic 
Levy variables and quantify the dependence of x and v- 
variables. 



II. GENERAL CONSIDERATIONS 

Although the initial part of our calculations does not 
go beyond fl2l |. we reproduce part of the considerations 
in a physical language for the sake of readability. The 
formal solution of Eq. ([]} is 



x{t) = F(t) 



G(t - t')£(t')dt' , 



(3) 



where G(t) is the Green's (response) function of the cor- 
responding process, and F(t) is a decaying function (a 
solution of the homogeneous equation under given initial 
conditions). The solution for v is given by 



v{t) = F v (t) + / G v {t - t')Z(t')dt' , 



(4) 



where G v (t) is the Green's function of the velocity pro- 
cess 



G v (t) = -G(t). 



(5) 



In a stationary situation, t — > oo, the F- functions in 
Eqs. ([3]) and (jU vanish. The Green's function of Eq. ([T]) 
can be easily found e.g. via the Laplace representation, 
and reads: 



_ ex P (- 7 V2) _ 



V" 2 - 7 2 /4 
for uj = \f~k > 7/2 (underdamped case 
G(t) =texp(-jt/2) 
for uj = 7/2 (critical case) and 



\/w 2 -7 2 /4i 



vV/4- w 2 t 



(6) 



(7) 



(8) 



for uj < 7/2 (overdamped case). Note that the functions 
G(t) vanish both for t = and for t — > 00 so that 



G(t) 



(If 



G(t) 



dt 



1 

1 



(/f 



G 2 W 



dt (9) 



= o G Wlo -0, 



i.e. G(t) and G«(t) are orthogonal on [0,oo). 

Let us now confine ourselves to the finite (but long 
enough) interval of integration and consider integral sums 
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FIG. 1. Joint probability density p [x, v) (top panel) and prod- 
uct of two marginal densities p x (x)p v (v) (bottom panel) for 
a = 1. The left panel corresponds to the underdamped case 
(k = 1,7 = 1) while the right one to the overdamped case 
(k — 1,7 = 4). Dashed lines in the top panel represent major 
axes described in the text. Simulation parameters: number 
of repetitions N = 10 7 , time step of integration At — 1CP 2 . 



corresponding to Eqs. ([3]) and ((4]) for a given time step 
At: 



N 



(10) 



and 



JV 



(ii) 



1=1 



with £j being independent identically distributed random 
variables corresponding to the integrals of the noise over 
independent At- intervals, and flj, hi being the values of 
the Green's functions evaluated at corresponding points. 

Let us evaluate the joint PDF of x and v 



p(x,v) 



(12) 



The characteristic function of this distribution is given 
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by the Fourier transform 

f(k,q) = [[ e lkx+lqv p(x,v)dxdv 



(13) 



N 



z=l J 
N 

~[f(ka l + qb i ), 



where f(k) = J e tk tp(£)d!; is the characteristic function 
of the distribution of £. 

In the case considered here, the noise corresponds to 
a formal derivative of a Levy random process L(t). The 
integrals of the noise over finite time intervals are in- 
dependent Levy stable random variables. The probabil- 
ity distributions corresponding to such random variables, 
stable infinitely divisible laws, are limiting distributions 
in the space of probability densities due to generalized 
central limit theorems 0, [lj|. The characteristic func- 
tion of £(Ai) is thus given by 



D ife«(At)\ 



exp 



-Ata a \k\ a (l-i/3sgnfctan 



7ta\ 



(14) 

where a £ (0,2] is the stability index, j3 £ [—1,1] is the 
asymmetry (skewness) parameter, while a (a > 0) is a 
scale parameter. Increments of the Levy motion are dis- 
tributed according to a a-stable density, which for a < 2 
has power law asymptotics, i.e. p(£) oc 1/(£) Q+1 - The 
Gaussian distribution is a special case of a stable den- 
sity with a = 2, which is the only one a-stable density 
possessing finite moments [l[. In our work here we con- 
centrate on the case of symmetric distributions, /3 = 0. 
The characteristic functions of symmetric Levy distribu- 
tions have the form 



f(k) = exp(-Ato"|fc| a ). 



Therefore 



'J f(kdi + qbi) = cxp 



i=i 



N 



-Ata a ^2\ka t + qbi 



or returning to the integral notation, 



f(k, q) = exp 



\kG(t) + qG v (t)\ a dt 



(15) 



(16) 



(17) 



The characteristic functions of the marginal distributions 
are obtained by putting q = (for ^-distribution) or 
k = (for w-distribution) . Note moreover, that the char- 
acteristic function of the distribution of any linear com- 
bination of x and v 



will read 



z = ax + bv 



f z (k) = f(ak,bk) 



(18) 
(19) 



and is the one of the symmetric Levy distribution, so that 
the distribution of x and v is a bona fide bivariate stable 
one. It is however not of elliptic type, for which case the 
characteristic function would be 



f(k,q) = cxp 



2;„2 



C xv kq 



wW\ a ^dt 



20) 



i.e. a function of a quadratic form in k and q, see [If 

The characteristic function, Eq. (fl~7)) does not decouple 
into the product of the functions of k and q only, i.e. the 
distributions of x and v are not independent (we shall 
say, they are associated) , except for the Gaussian (a = 2) 
case. In this case 



\kG{t) + qG v {t)\ z dt = k A G 2 (t)dt (21) 

Jo 



/•OO 

+ q 2 G 2 v (t)dt 
JO 

/>oo 

+ 2kq / G(t)G v (t)dt, 
Jo 



and the last integral vanishes due to orthogonality of G(t) 
and G v (t). In all other cases such decoupling does not 
take place. The Gaussian case is also the only case when 
the distribution is elliptic. 

In order to demonstrate lack of independence of x and 
v for the harmonic oscillator under Levy noise we have 
estimated from numerical simulations of the process the 
joint distribution p(x, v) and compared it with product of 
two marginal densities p x (x)p v (v), see Fig. [TJ presented 
both for the underdamped (left panel) and for the over- 
damped (right panel) case under Cauchy (a = 1) noise 
[13] ■ The contours of the corresponding joint probabil- 
ity densities are central symmetric, and the positions of 
their main axes, along which the width is the largest or 
the smallest, can be obtained numerically by finding the 
extrema of 



/* OO 

j(9) = / \G(t) cos + <?„(*) sin ( 
Jo 



'dt. 



(22) 



The corresponding axes are at the angles of approxi- 
mately —43° (longer axis) and 47° (shorter one) for the 
underdamped case, and of approximately —15° (longer 
axis) and 75° (shorter one) for the overdamped case. 
The corresponding angles are close (but not equal to) the 
positions of extrema in the phase plot, due to different 
weighting of the variables, compare Figs. [1] and [3] 



III. MEASURES OF DEPENDENCE 

The case of dependent random variables lacking mo- 
ments poses a problem of quantifying the strength of their 
dependence. In the "normal" Gaussian case, and in every 
case when the second moments of the variables are finite, 
the coefficient of covariance gives us a standard measure 
of such dependence. For variables with zero mean 



cov(x,y) = (xy), 



(23) 
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and for variables whose mean values are nonzero the same 
is defined for the centered variables and 
y c = y — y. In what follows we consider only variables 
centered at zero, and omit the subscript c. To character- 
ize the relative strength of the dependence one can intro- 
duce the correlation coefficient corr(x, y) by normalizing 
the covariance over the dispersions of the corresponding 
centered variables. Alternatively, one can consider the 
covariance of variables x and y normalized with respect 
to their dispersions a x = (i 2 ) 1 / 2 and a y = (y 2 ) 1 ^ 2 , i.e. 
x = x/a x , y = y/a y : 

, . cav(x,y) ,~ ~\ 

corr(a;, y) = = cov(x,y). (24) 

<J X (7y 

The correlation coefficient of two random variables 
changes in the interval [—1, 1] and it is unity if the vari- 
able y is a copy of x up to arbitrary rescaling, and to 
minus unity if y is a copy of x taken with the oppo- 
site sign, again up to the arbitrary change of scale. The 
correlation coefficient of independent variables vanishes. 
Since the definition of the correlation coefficient involves 
the second moment of the corresponding random vari- 
ables, it cannot be used immediately for the ones lacking 
such a moment. Here several generalizations have been 
proposed, all seeming to be deficient for our purpose, 
consequently we have to extend their definitions. 

A relatively common measure of dependence of Levy 
variables is the codifference fl8| . To understand its na- 
ture it is enough first to consider the Gaussian case and 
to turn to the characteristic function of the corresponding 
bivariatc distribution of x and y 

f(ki,k 2 )=[ J dxdye^ x+ ^y p (x,y) (25) 
= cxp [— u 2 fc 2 — cryk 2 — 2cov(cc, y)k\k2\ ■ 

To assess the covariance one can then consider the dif- 
ference variable d = x — y whose characteristic function 
is 

f d (k) = f(k,-k) (26) 
= cxp [- (a 2 x + a 2 y - 2cov(a;, y)) k 2 ] , 

from which 

cov(z, y) = X - [In f d (k)\ k=1 + o\ + a 2 x \ , (27) 

or in slightly different notation 

cov(a:, y) = -codiff (x, y) (28) 

with codiff (x, y) being the codifference of the variables x 
and y defined through 

codiff(z )2 /) =lnf(k,-k)\ k=1 (29) 
+ ln/(fc,0)| fc=1+ ln/(0,fc)| fc=1 . 

The definition of the codifference as a measure of depen- 
dence via characteristic functions does not rely on the 



existence of any moments and is universal, i.e. it can be 
used in the Levy case as well. We note that equivalently 
one can use the sum of the two variables x and y, and to 
define the "cosum" as 

cosum(a;, y) = — In f(k, k)\ k=1 (30) 

+ ]n/(M)| fc=1 +]n/(0 I ft)| fc=1 

The codifference (or cosum) as a measure of dependence 
has two important drawbacks, the first one quite specific 
to their definitions, and the second one inherited from 
the the definition of covariance itself. 

First, defining the corresponding parameters via char- 
acteristic functions, being the means of strongly oscil- 
lating exponentials implies numerically generating very 
large data samples. Therefore, for the sake of practica- 
bility, the measure of dependence has to rely on some 
robust and easily accessible statistics. 

Second, using the codifference as a measure of depen- 
dence of the observablcs of different dimension (coordi- 
nate and velocity) makes no sense physically, since the 
change in units of measurement changes the measure of 
dependence. The codifference measure is often used in 
the analysis of time series, since the dimensions of all 
terms are the same, the problem of units therefore do 
not stand. The same is true also for the normalized cod- 
ifference (liij ). which is to no extent a simple analogue of 
the correlation coefficient. 

Let us first address the problem of statistics. For a 
multivariate Levy distributions, as the one appearing in 
the previous Section, both the distributions of the the 
difference d = x — y, pd(d), and of the sum s = x + y, 
p s (s), of x and y (in our initial problem the variable y 
corresponds to the velocity v of the particle) are univari- 
ate Levy distributions. The widths of these distributions 
are exactly the prefactors of |fc|" in the corresponding 
characteristic functions 

< = |ln/(fe,-fe)| fc=1 | (31) 

and 

< = |ln/(fe,fe)| fc=1 |. (32) 
Analogously w x and w y are equal to 

< = |ln/(fc,0)| fe=1 | (33) 

and 

< = |ln/(0,fc)| fc=1 |. (34) 

Therefore 

codiff (a;, y) = w°[ - w% - w^, (35) 

and 

cosum(.T, y) = — w x — Wy . (36) 

Contrary to the initial definition via characteristic 
functions, the ones over the widths (scaling parameters) 
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of the distributions can be easily assessed numerically 
since they are connected with the robust characteristics 
like the interquartile distances of the corresponding dis- 
tributions or with the heights of their probability densi- 
ties at zero 



»(0) - 1 



—w 01 \k\' 



dk = 



r(i/«) 



TTW 



(37) 



which are easier to read out from the numerically ob- 
tained histograms. This "height measure" is repeatedly 
used in our work. 

To circumvent the second problem, we propose the fol- 
lowing measure of dependence D corresponding to the 
cosum of the normalized variables. Both distributions, 
the one of x and the one of v are first normalized by in- 
troducing x = x/w x and v = v/w y (both possessing Levy 
distributions of unit width) . Then we compare the width 
of the distribution of C = x + v (raised to the power of 
a) with 2, the sum of width of the distributions of x and 
v: 



D = 



(38) 



and call this measure of dependence the width excess (the 
letter D stands here for "Dependence"). The measure D, 
the cosum of rescaled variables, changes in the interval 
[—1,1]. It is zero for independent variables, equal to 1 
if the variable v is a copy of x up to arbitrary rescaling 
(i.e. up to the change of units) and is equal to —1 if the 
variable v is a copy of —x, up to rescaling. It is clear 
that D is insensitive to the units in which x and v are 
measured. 

For a bivariate Gaussian density D is exactly the cor- 
relation coefficient 



D = 



= corr(i, y). 



(39) 



For the elliptic Levy process it is connected with the 
Press' association parameter [l6| . which however can 
hardly be generalized to nonelliptic cases and vanishes 
for isotropic distributions, the ones with characteristic 
function 11811 



/(fc, q) = exp 



-a a (k 2 +q 2 )^ 2 



(40) 



although the corresponding x and v values are not inde- 
pendent, see Eqs. (J35|) and (T4"4"|). 

The characteristic function of the distribution of £ 

£ = x + v = x/w x + v/w v 



IS 



ft( k ) = f{k/w x ,k/w v ) 
Consequently the width excess is 

G(t) , G v (t) ' 



D= 1 - 
2 



h 



dt-l 



(41) 



with 



and 



h = 



h = 



\G{t)\ a dt 



\G v (t)\ a dt 



l/a 



1/a 



(42) 



(43) 
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FIG. 2. Width excess D, see Eq. (|41|) . for various a. Lines 
represent theoretical formulas for underdamped (fc = 1,7 = 
1) and overdamped (k = 1,7 = 4) cases, points correspond 
to simulation results. The dashed line presents values of the 
width excess D for isotropic case, see Eqs. (|40|l and 
Other simulation parameters as in Fig. [1] 



The width excess (|41j) vanishes for a ~ 2 only, in- 
dicating the independence of the position and velocity 
variables. Contrary to the a = 2 case, for a < 2 it takes 
negative values, indicating anti-association, i.e. depen- 
dence between the position and the velocity. The results 
of simulations and calculations of the width deficit for the 
underdamped and overdamped cases are given in Fig. [21 
together with the result for an isotropic case, for which 



(44) 



D = 2"/ 2 - 1 - 1 



corresponding to all elliptic distributions, see Eqs. (|38|) 
and (|40|) . Figure [2] demonstrates perfect agreement be- 
tween theoretical results given by formula (|4ip (lines) and 
numerical data (points). Negative values of the width 
excess D indicate strong negative association, which is 
stronger for the overdamped case (k = 1,7 = 4) than for 
the underdamped case (k = 1,7= 1). 



IV. PHASE AND AMPLITUDE 
DISTRIBUTIONS 

The strong anti-association of velocity and coordinate 
is reflected in the distribution of the phase angle 



arctan 



arctan 



v/w v 



X/Wa 



(45) 
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FIG. 3. Phase distribution p{(j>) in the underdamped (k — 
1,7 = 1- top panel) and overdamped (k = 1, 7 = 4 - bottom 
panel) cases for different values of a. The full line shows the 
corresponding distribution for independent Cauchy variables, 
see Eq. (|46l) . Simulation parameters as in Fig. Q] 



This distribution is homogeneous over [— 7r/2,7r/2) for 
all elliptic bivariate Levy laws (e.g. in the Gaussian 
case) and has a marked form with peaks (divergences) 
at <f> — 0, ±7r/2 for independent random Levy variables 
(corresponding to the star-like form of the product of two 
Levy distributions, see Fig. [IJ. The actual forms of the 
phase PDFs are shown in Fig. [3] for the underdamped 
(top panel) and overdamped (bottom panel) cases. It 
is interesting that the peaks are present at "nontrivial" 
angles, the ones different from and ±7r/2. 

For a < 2, strong deviations both from the uniform 
distribution and from the product distribution are well 
visible. For a = 1 the phase distribution for independent 
x- and v-variables is given by 



p{4>) = 

see Fig. El 

The amplitude 



7T 2 cos2</> 



In 



1 + cos 2<j> 



1 — cos 2(f> 



(46) 



A = \Jx 2 + v 2 = ^f(x/w x ) 2 + {v/w v ) 2 (47) 
is asymptotically distributed according to a power-law 



FIG. 4. Complementary cumulative distribution (1 — F(A)) 
of the amplitude A — V x' 2 + v 2 = J (x/w x ) 2 + (v/w v ) 2 in 
the underdamped (k = 1,7 = 1 - top panel) and overdamped 
(k — 1,7 = 4 - bottom panel) cases for different values of 
a. Full lines show the corresponding asymptotic power-law 
decays of complementary cumulative distributions, i.e. A~ a . 
Simulation parameters as in Fig. [T] 



distribution 



p(A) oc 



A a+1 ' 



(48) 



consequently the complementary cumulative distribution 
has also power-law tails which are characterized by the 
exponent a, see Fig. 21 as it follows from the results of 
Ref. [ll. 



V. SUMMARY 

In situations far from equilibrium interactions of the 
harmonic oscillator with the environment can often be 
described by the white Levy noise. Presence of Levy 
noises introduces dependence (association) between ve- 
locity and position, which vanishes only in the Gaus- 
sian limit of a = 2. The presence of association be- 
tween position and velocity is manifested by the nontriv- 
ial joint distribution of p(x, v) and nontrivial phase dis- 
tribution p(4>). Our main finding corresponds to strong 
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anti-association between the position and the velocity of 
a Levy-driven oscillator, a property which might be of 
high importance for first passage properties of such a 
process. For example, in the Kramers problem the de- 
pendence of position and velocity might be responsible 
for the breakdown of the transition-state description. A 
direct consequence of (anti-)association between velocity 



and position is the breakdown of equipartition theorem 
of equilibrium statistical physics for the harmonic Levy 
oscillator. The recorded findings, together with results of 
earlier investigations 0, [2(| [21[ clearly show that prop- 
erties of out-of-equilibrium (Levy noise driven) systems 
are very different from their equilibrium (Gaussian noise 
driven) counterparts. 



[1] A. Janicki and A. Weron, Simulation and chaotic behav- 
ior of a- stable stochastic processes (Marcel Dekker, New 
York, 1994); A. V. Chechkin, V. Y. Gonchar, J. Klafter, 
and R. Metzler, Adv. Chem. Phys. 133, 439 (2006); A. A. 
Dubkov, B. Spagnolo, and V. V. Uchaikin, Int. J. Bifur- 
cation Chaos. Appl. Sci. Eng. 18, 2649 (2008). 

[2] M. F. Shlesinger, G. M. Zaslavski, and J. Klafter, 
Nature (London) 363, 31 (1993); J. Klafter, M. F. 
Shlesinger, and G. Zumofen, Phys. Today 49, 33 (1996); 
T. H. Solomon, E. R. Weeks, and H. L. Swinney, 
|Phys. Rev. Lett] 71, 3975 (1993); D. del Castillo- 
Negrete, Phys. Fluids 10, 576 (1998). 

[3] A. V. Chechkin, V. Y. Gonchar, and M. Szydlowski, 
Phys. Plasmas 9, 78 (2002); D. del Castillo-Negr ete, 
B. A. Carreras, and V. E. Lynch, Phys. Re v. Lett] 94, 
065003 (2005). 

[4] H. Katori, S. Schlipf, and H. Walther, |Phys. Rev. Lett] 
79, 2221 (1997). 

[5] C.-K. Peng, J. Mietus, J. M. Hausdorff, S. Havl in, H. E. 
Stanley, and A. L. Goldberger, |Phys. Rev. Lett.| 70, 1343 
(1993). 

[6] R. Segev, M. Benveniste, E. Hulata, N. Cohen, 

A. Palevski, E. Kapon, Y. Shapira, and E. Ben-Jacob, 

|Phys. Rev. LetE| 88, 118102 (2002). 
[7] M. A. Lomholt, T. Ambjornsson, and R. Metzler, Phys. 

Rev. Lett. 95, 260603 (2005). 
[8] G. M. Viswanathan, V. Afanasyev, S. V. Buldyrev, E. J. 

Murphy, P. A. Prince, and H. E. Stanley, Nature (Lon- 



don) 381, 413 (1996). 
[9] P. D. Ditlevsen, Geophys. Res. Lett. 26, 1441 (1999). 

[10] R. N. Mantegna and H. E. Stanley, An introduction 
to econophysics. Correlations and complexity in finance 
(Cambridge University Press, Cambridge, 2000). 

[11] D. Brockmann, L. Hufnagel, and T. Geisel, Nature (Lon- 
don) 439, 462 (2006). 

[12] B. J. West and V. Seshadri, Physica A 113, 203 (1982). 

[13] J. L. Doob, Ann. of Math. 43, 351 (1942). 

[14] S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943). 

[15] B. V. Gnedenko and A. N. Kolmogorov, Limit dis- 
tributions for sums of independent random variables 
(Addison- Wesley, Reading, MA, 1968). 

[16] S. J. Press, J. Multivariate Anal 2, 444 (1972). 

[17] P. Garbaczewski and R. Olkiewicz, J. Math. Phys. 41, 
6843 (2000). 

[18] G. Samorodnitsky and M. S. Taqqu, Stable non-Gaussian 
random processes: Stochastic models with infinite vari- 
ance (Chapman and Hall, New York, 1994). 

[19] D. Rosadi, Stat, and Prob. Letters 77, 1490 (2007). 

[20] A. V. Chechkin, J. Klafter, V. Y. Gonchar, R. Metzler, 
and L. V. Tanatarov, Phys. Rev. E 67, 010102(R) (2003). 

[21] A. V. Chechkin, V. Y. Gonchar, J. Klafter, R. Metzler, 
and L. V. Tanatarov, J. Stat. Phys. 115, 1505 (2004). 



